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We demonstrate that the quasi-equilibrium state in self-gravitating JV-body system after cold col¬ 
lapse are uniquely characterized by the local virial relation using numerical simulations. Conversely 
assuming the constant local virial ratio and Jeans equation for spherically steady state system, we 
investigate the full solution space of the problem under the constant anisotropy parameter and ob¬ 
tain some relevant solutions. Especially, the local virial relation always provides a solution which 
has a power law density profile in both the asymptotic regions r — > 0 and oo. This type of solu¬ 
tions observed commonly in many numerical simulations. Only the anisotropic velocity dispersion 
controls this asymptotic behavior of density profile. 

PACS numbers: 98.10.+z,05.10.-a 

I. INTRODUCTION 

and clusters of galaxies are typical self¬ 
systems (SGS) in the universe. Extracting 
the essence of them, study on the quasi-equilibrium state 
of SGS is an important clue for understanding the basic 
properties of structures in the universe. In the process 
of the structure formation of SGS, a cold collapse and 
a cluster-pair collision would be the most fundamental 
processes. Therefore in this paper, we would like to focus 
on the characteristic properties of the quasi-equilibrium 
state formed by such basic processes. 

For many numerical simulations, after a cold collapse, 
system settles down to a quasi-equilibrium state which 
is well described by the Vlasov equation for a spheri¬ 
cally symmetric steady state system. In such a quasi¬ 
equilibrium state, the spherically averaged densitypro- 
file is found to be p oc r -4 at the outer regionjlj, 2 
and the velocity dispersion is isotropic near the cen¬ 
tral region and becomes radially anisotropic at the outer 
regionj3|. Recently we found two characteristic proper¬ 
ties which appearer in the quasi-equilibrium state after 
a cold collapse 0 Of]. One is the linear temperature-mass 
relation which yields a characteristic non-Gaussian veloc¬ 
ity distribution. The other is the local virial (LV) relation 
which is robust against various initial condition such as 
the density profile and the virial ratio. 

The Plummer model is one of the examples which 
obey the LV relation. N.W.Evans and J.An discuss 
the anisotropic distribution function where the LV re¬ 
lation is satisfied ( they call it hypervirial ) and obtain 
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the analytical solution which has a constant anisotropy 
parameter |(|. In this paper, we would like to explore 
the relevance of LV relation in dynamics of SGS. There¬ 
fore, we assume the general LV relation where the virial 
ratio is constant and investigate the role of LV relation 
in a spherically symmetric steady state system. More¬ 
over, we compare the solutions of Jeans equation under 
the LV relation and the results of numerical simulations. 
In numerical simulations, we use a leap-frog symplectic 
integrator on GRAPE-5, a special-purpose computer de¬ 
signed to accelerate N-body simulations^)- 

In section HU we first show the validity of LV relation 
with the numerical N-body simulation. Then, assuming 
the constant virial ratio, we investigate the full solution 
space of Jeans equation under the constant anisotropy 
parameter in section eh We show the special properties 
of the LV relation and obtain analytically the physically 
relevant solutions. In section El we compare these solu¬ 
tions with the results of numerical simulation and discuss 
the validity and the consistency of the LV relation. The 
last section m is devoted to the discussions and further 
developments of the present work. 

II. LOCAL VIRIAL RELATION 

In a steady state, the gravitationally bound system 
settles down to a virialized state satisfying the condition 

W + 2K = 0, (1) 

where W and K are, respectively, the time-averaged po¬ 
tential energy and kinetic energy of the whole bound sys¬ 
tem. This is a well-known global relation that holds for 
the entire system. 

Here we define the LV relation at each position r as 
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FIG. 1: The LV relation for some numerical simulations ob¬ 
tained by a typical cold collapse simulation; spherical col¬ 
lapse (SC) and cluster-pair collisions (CC). In the case of SC, 
5000 particles are distributed with a power law density profile 
(p oc r~ a ) within a sphere of radius R and the initial virial 
ratio (bt) is set to be small. In the case of CC, each cluster 
has the equal number of particles (2500) and all particles are 
homogeneously distributed within a sphere of radius R and 
is set to be virialized initially. The initial separation of the 
pair is 6 R along the x axis. In all of the simulations, soft¬ 
ening length e = 2 ~ 8 R is introduced to reduce the numerical 
error caused by close encounters. The LV ratio b is plotted 
as a function of M r /Mt.ot, where Mtot is the total mass of the 
system. The virial ratios at M r are time averaged. 


where a and $ are the velocity dispersion ( a 2 = of + 
0 g + aV) and the potential energy, respectively. Hence, 
we define the LV ratio 


b(r) := -2o- 2 (r)/3>(r), (3) 


in order to measure to what extent the LV relation 
holds. Actually, this LV relation appearers in the quasi¬ 
equilibrium state obtained by some numerical simula¬ 
tions iS For a cold collapse and cluster-pair collision 
from a initially smooth density profile, the LV ratio is 
shown in FiglU The value b(r) takes almost unity in the 
inner region and slightly reduces in outer region for all 
of the simulations. On the other hand, for the case of 
SC (bt = 0, a = 2) where the initial virial ratio b t = 0 and 
the initial density profile p oc r~ 2 , the deviation at the 
central region is significant. 

For a wide class of collapses from various initial condi¬ 
tions such as the change of the virial ratio and the density 
profile, the LV relation are almost satisfied. In the next 
section, we investigate the role of the LV relation in a 
steady state of the SGS. 


III. LOCAL VIRIAL RELATION AND JEANS 
EQUATION 


To begin with, we assume the general LV relation, 


2 q - 2 (r) , 

Hr) 


(4) 


where the relative potential </>(r ) := — $(r) and b is a 
constant. When 6 = 1, LV relation is exactly satisfied. 
The Jeans equation for a spherically symmetric steady 
state system is 


d(pa^) 2/3 2 _ d<t> 

—1 - 1 - P a r ~ P~T~ 

dr r dr 


where the anisotropy parameter /3 is defined as 


(3~1- 
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(5) 

( 6 ) 


We study the density profile which satisfies the general 
LV relation and Jeans equation for a spherically sym¬ 
metric steady state system ®. We consider the sim¬ 
plest case, where (3 is a constant. Following W. Dehnen 
and D.E. McLaughlin ;;S], we investigate the full solution 
space of the problem. 

Eliminating a 2 from Eq. 0 by using Eq. ®, the rela¬ 
tive potential becomes 


c/> = cf)ox 2 P/ a y V“, (7) 

where x := r/ro, y := p/po, a := (6 — 4/3— b)/b, and <f >o := 
</>(ro). Substituting Eq.0 into the Poisson equation, the 
density profile satisfies the following equation: 

i - -(7 - 7a)(7 - 7b) = aKX 2 ~ 20/a y 1 ^ 1/a , (8) 

a 

where a prime denotes the differentiation with respect to 
lnx and 7 := — dlny/dlnx and k := ATTGr^po/ipo, and 

7 a ■= 2/3, 7 b := 2/3 + a. (9) 

We will restrict our consideration to the case with /3 < 1 
and 0 < b < 2 . 

First, the equation 0 has a scale invariant solution 
y = x~ lx with 


7i 


2(a ~ (3) 

a — 1 


( 10 ) 


for a ^ 1. In the case for a = 1, Eq.0 has no scale 
invariant solution except for (/3, b) = (1,1). For (/3,6) = 
(1,1), the r.h.s of Eq.0 becomes always constant (k) and 
we can easily integrate the Eq.0 and the total mass is 
infinite ( see Appendix fXll. Hereafter we consider the case 
for a/1 (we discuss the case for a = 1 in Appendix 0J. 
In this case, the structure of the phase space ( 7 /, 7 ) is 
classified into three cases. 
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TABLE I: The asymptotic behavior of 7 in Eq. ©. 
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(B) 7a < 71 < 76 for 

3 — 26 

(3 < 1(0 < 6 < 1) or (3 < ———(1 < 6 < 2), 

(C) 7a < 7fe < 7i for 

3-26 3-6 
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truncation. This second family is also unphysical because 
this family has negative 7 b Finally, the solution which 
is represented by a sold line in Fig© is realistic, because 
the density profile behaves p oc r~ la as r —> 0 and has 
an outer truncation beyond a finite radius. 



For each case, the asymptotic behavior of 7 is sum¬ 
marized in Table [I] Note that the solution with the 
power law behavior in both of the asymptotic regions 
(r —> 0,oo) exists only in the case B. Since n > 0 and 
a > 0, non-negativity of the r.h.s of Eq. © leads to the 
following condition: 

7'>-(7-7a)(7-7b)- (11) 

a 

From the condition CD. both of the case A and case C 
have no scale invariant solution ( 71 ). 

Now investigate the solutions of Eq. © in detail. Dif¬ 
ferentiating Eq. © with respect to lnx, we get the fol¬ 
lowing closed second order differential equation for 7 , 


7 


// 


3 — a . 

-7 

a 


7 - o-(7a + 7b - (a - 1 ) 71 ) 

O OL 
OL — 1 

= —2 - (7-7a)(7-7b)(7-7i)- (I 2 ) 
or 


Using Eq. o, we study the flow of solutions in the ( 7 ', 7 ) 
phase space for the above three cases. Note that the 
solution with 7 ' < 0 is unphysical because this solution 
is unstable. 

In the case A, the flow of the solutions in Eq. Cl is 
shown in Fig© All solutions have an outer truncation 
beyond a finite radius (7 —> 00 ). The behavior of the 
solutions is classified into two families by the solution 
which starts from the fixed point 7 = and is repre¬ 
sented by a solid line in Fig© Solutions which exist in 
upper region of this solid line, are that the density pro¬ 
file has an inner hole and an outer truncation. These are 
unphysical solutions because the solutions with an inner 
hole are unstable. On the other hand, solutions which 
exist in lower region of this solid line, are that the den¬ 
sity profile behaves p oc r~ lb as r —> 0 and has an outer 


FIG. 2: The flow diagram in the (' 7 ', 7 ) phase space for case 
A. The parameters are set to (/3, b) = (7/8, 3/2). A dashed 
line represents the boundary defined by Eo. ll It . Two filled 
circles denote the fixed points 7 a (left) and 76 (right), respec¬ 
tively. 


In the case C, the behavior of the ( 7 ', 7 ) phase space 
is shown in Fig© and is same as one in the case A. The 
realistic solutions are represented by a sold line where 
the density profile behaves p oc r~ la as r —> 0 and has 
an outer truncation. 
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FIG. 3: The same as Fig© but for case C. The parameters 
are set to (/3, b) = (1/2, 3/2). 
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Figs©© and©show the behavior of the ( 7 7 , 7 ) phase 
space in the case B. In this case, the behavior of the 
solutions in Eq. d depends on whether 6 is less than, 
greater than, or equal to a critical value 6 cr q = 1 , where 
71 = (q a + 7b)/2 and the LV relation is satisfied. 

For 6 < 6 cr p, the flow of the solutions is shown in 
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FIG. 4: The same as FigEl but for case B (0 < b < 1). The 
parameters are set to (/3, b) = (3/4, 3/4). Three filled circles 
denote the fixed points -y a (left), 71 (middle), and 75 (right), 
respectively. 


FIG. 6 : The same as FigEl but for case B (b = 1). The pa¬ 
rameters are set to (/3, b) = (1/2,1). The solid line connected 
between 7 a and 7 b denotes the critical solution. 
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FIG. 5: The same as FigEl but for case B (1 < b < 2 ). The 
parameters set are to (/3, 6 ) = (5/7, 9/8). 


FigEl In this case, all solutions except one are unphysical 
because the density profile has an inner hole. The only 
exceptional solution is represented by a solid line which 
starts from 7 = j a and ends at 7 = 71 in FigEl The 
density profile of this solution behaves as p oc r _7a in the 
limit of v —► 0 and undergoes damped oscillation around 
p oc r -71 in the limit of r — > 00. 

For b > b CT the flow of the solutions is shown in FigEl 
The behavior of the solutions in this case corresponds 
to a mirror image of one in the case b < 6 cr p. The 
behavior of the solutions is classified into two families by 
the solution which starts from the fixed point 7 = y a and 
is represented by a solid line in FigEl Solutions which 
exist in upper region of this solid line, are unphysical 
because the solutions have an inner hole. On the other 
hand, solutions which exist in lower region of this solid 
line, go around p oc r -71 and have an outer truncation 
beyond a finite radius or end at 7 = 7& (p oc r ~ lb ). This 


second family with positive 7 ' is physical. Similarly, the 
solution which starts from 7 = 7 o and is represented by 
a sold line in FigEl is also realistic, because the density 
profile behaves p oc r _7a as r —> 0 and has an outer 
truncation beyond a finite radius. 

FigEl shows the how of the solution for b = ^ 7 .^. In 
this case, there is a following first integral, 


I< = 


, 2(1-13). v J 2(1H3) 

7 + 5 — 4/3 W ~ 7a ^ 7 ~ lb > 

1 


7' - 5 _ 4 ^ (7 - 7a) (7 - 7b) 

Especially, K = 0 leads to the Eq. m and 

2 ( 1 -/?), 


7 = 


5-4/3 


■(7 — 7a) (7 — 7b)- 


(13) 


(14) 


The above equation o shows the critical solution, 
which starts from 7 = 7 a and ends at 7 = 7 b and is 
represented by a solid line in FigEl Solutions which ex¬ 
ist in upper region of this critical solution, have an inner 
hole. Others which exit in lower region of the critical so¬ 
lution, eternally oscillate around the fixed point 7 = 7 i- 
From these behaviors, all solutions are unphysical except 
the critical solution. The only critical solution is physical 
and has a special characteristic that it shows power law 
behavior in the both limiting region r —> 0 (p oc r _7a ) 
and 00 (p oc r~ lb ). 

It is easy to obtain an analytical form of the critical 
solution. From Eq. we get 


7 = 


7 a + 7 bX 


2 ( 1 — 0 ) 


l +x 2(l-P) 


(15) 


where we choose an integral constant as 7 ( 1 ) = 71 . Inte¬ 
grating again Eq. ©> we obtain 


P = ^jM tot x 2 - s ( l+x s )~^ +1 / s \ (16) 
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M = M tot { l+x~ s )-( 1+1 ' s \ (17) 

<t> = ——^(1 + x s )~ 1/f/l , (18) 

ct2 = GM» I/., (ig) 

zr 0 

^ ») 

where s := 2(1—/?), M := f* An^pctu, and M to t is a total 
mass (47rrpp 0 )/(l+s). This critical solution was found by 
U.-I.K. Veltmann 0 and N.W.Evans and J.An obtained 
by solving Jeans equation assuming the LV relation ||. 

The role of the LV ratio is similar to the scale invari¬ 
ant phase space density which is observed in the cosmo¬ 
logical simulations based on the cold dark matter sce¬ 
nario. Assuming the scale invariant phase space density 
(<t 3 //9 oc r~ a ) and solving Jeans equation for spherically 
steady state systemQ, the critical value a cr = 35/18 ex¬ 
ists and the flow diagram in the phase space around the 
critical value is similar to one in the case with the LV ra¬ 
tio. The critical solution with the anisotropy parameter 
(3 = 7/8 shows the scale invariant phase space density 
(a = 7/4 < a cr ). 

IV. COMPARISON WITH SIMULATION 

In the previous section, we showed that the LV relation 
is satisfied quite well for cold collapse simulations except 
for the case with a highly concentrated initial matter dis¬ 
tribution. We can classify the N-body simulations into 
two main classes from the viewpoints of the behavior of 
anisotropy parameter /? for the bound state after a cold 
collapse. First is the ones with the initial homogeneous 
sphere where the particles are distributed homogeneously 
in the spherical region, which we call class I. In these ini¬ 
tial conditions, large amount of particles reach the mass- 
center at the same time, which causes the high density 
central region with the isotropic velocity dispersion after 
the collapse. Actually we got the results that (3 vanishes 
in the region within the central 30 percent cumulative 
mass regardless of the initial virial ratio in such cases 
(Fig0. On the other hand, the bound state starting 
from the other initial conditions show that (3 is positive 
and monotonically increases against the cumulative mass 
(Fig0. Hence we classify the simulations with these ini¬ 
tial conditions into the second class and call it class II. 

Here we compare the results of these N-body simula¬ 
tions with those of the critical solutions in the previous 
section from the viewpoints of mass density, temperature- 
mass (TM) relation, and phase space density. In cosmo¬ 
logical simulations, it is well known that the phase space 
density has a remarkable character, i.e., it has a single 
exponent with radius in the wide range of the viriallized 
region, despite the fact that the density has two different 
exponents El However, it is not obvious that it be¬ 
haves in the same way for the simulations with different 



FIG. 7: The distribution of the anisotropy parameter (3(r) for 
the simulations of class I. A plot with error-bar is the result 
of the numerical simulations SC with ( bt,a) = (0,0) (gray 
triangle) and ( bt,a ) = (0.5,0) (black diamond) in FigO A 
solid horizontal line represents the critical solution m-m 
with constant /3. Each line represents the critical solution 
with (3 = 0 (solid line), 0.25 (dot-dashed line), 0.5 (dashed 
line), and 0.75 (dotted line), respectively. 



FIG. 8: Same as Fig0 but for the simulations of class II. 
A plot is the result of the numerical simulations SC with 
( bt,a ) = (0,1) (black diamond), (0,2) (gray triangle), and 
CC (cross) in FigQ Each line represents the critical solution 
with (3 = 0 (solid line), 0.25 (dot-dashed line), 0.5 (dashed 
line), and 0.75 (dotted line), respectively. 


initial conditions. Hence it seems meaningful to examine 
the behavior of this quantity in our cold collapse simula¬ 
tions. 

First we compare the numerical results in the simula¬ 
tions of class I with the critical solutions with constant /?. 
The densities fit well to the distribution of the Plummer’s 
solution in the inner region as is shown in our previous 
paper (FigEJ. They, however, deviate from the Plum¬ 
mer’s distribution in the outer region where the velocity 
anisotropy is developed. We observe that the density in 
this region fits well rather to that of the critical solu¬ 
tion with positive constant (3. Both the TM relation and 
the phase space density have the same characters, i.e., 
both of them fit well to those of the Plummer’s solution 
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FIG. 9: A log-log plot for a density profile for class I simu¬ 
lations with the unit of rh = Mtot = G = 1 , where rh is the 
half-mass radius of the bound system. A plot is the result of 
the numerical simulations SC with ( bt , a) = (0, 0) (gray trian¬ 
gle) and ( bt,a ) = (0.5,0) (black diamond) in Fig0 Each line 
represents the critical solution with (3 = 0 (solid line), 0.25 
(dot-dashed line), 0.5 (dashed line), and 0.75 (dotted line), 
respectively. 



FIG. 10: Phase space density for class I simulations with the 
unit of ru = Mtot = G = 1. A plot is the result of the 
numerical simulations SC with ( bt,a ) = (0,0) (gray triangle) 
and ( bt,a ) = (0.5,0) (black diamond) in FigQ Each line 
represents the critical solution with (3 = 0 (solid line), 0.25 
(dot-dashed line), 0.5 (dashed line), and 0.75 (dotted line), 
respectively. 


in the inner region and to those of the positive constant 
/3 critical solution in the outer region (FigsITUI and EEU>. 
We can also derive the analytical solution which fits well 
to all of these physical quantities in the full region by 
connecting the inner Plummer’s solution with the outer 
critical solution with positive constant (3 ( AnnendixlBlor 
reference @). The phase space density in these simula¬ 
tions becomes flat in the central part as well as the mass 
density. In the critical solution, the central part is de¬ 
scribed by the fixed point 7 = 7 a , for which the exponent 
of the relative potential vanishes. Hence the exponent of 
the phase space density inevitably takes the same value 



FIG. 11: Temperature-Mass relations for class I simulations 
with the unit of rh = Mtot = G = 1. A plot with error-bar 
is the result of the numerical simulations SC with (bt,a) = 
(0,0) (gray triangle) and ( bt,a ) = (0.5,0) (black diamond) in 
FiglU Each line represents the critical solution with (3 = 0 
(solid line), 0.25 (dot-dashed line), 0.5 (dashed line), and 0.75 
(dotted line), respectively. 



FIG. 12: Same as FigEl but for class II simulations. A plot 
is the result of the numerical simulations SC with ( bt,a) = 
(0,1) (black diamond), (0,2) (gray triangle), and CC (cross) 
in Fig0 


as that of the mass-density at the center, as long as the 
phase space distribution follows the critical solution with 
constant (3. 

For the simulations of class II, the anisotropy parame¬ 
ter (3 behaves quite differently. In those cases the velocity 
dispersion is anisotropic even in the inner region. Actu¬ 
ally, the mass densities for those simulations are neither 
cuspy nor flat in the central region (FigE2J). Although 
the anisotropy parameter [3 is not constant but is mono- 
tonically increasing outward, the density profiles fit well 
to that of the critical solutions with /3=0.5 to 0.75 in full 
of the bound region. Both the phase space densities and 
TM relations for these simulations are also fitted to these 
critical solutions quite well, except for the case with the 
initial high density simulation with ( 6 t , a) = ( 0 , 2 ). 

In this last case with ( 6 t ,a) = (0,2), the phase space 
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FIG. 13: Same as Fip ll 01 but for class II simulations. A plot 
is the result of the numerical simulations SC with (bt,a) = 
(0,1) (black diamond), (0,2) (gray triangle), and CC (cross) 
in Fig[l| 



FIG. 14: Same as FigEl but for class II simulations. A plot 
with error-bar is the result of the numerical simulations SC 
with ( bt,a ) = ( 0 , 1 ) (black diamond), ( 0 , 2 ) (gray triangle), 
and CC (cross) in FigQ 


density becomes cuspy against the behavior of the mass 
density (Fig lldl) and the temperature falls off at the cen¬ 
ter against the behavior of any critical solutions (Fig IMfl . 
These results are certainly correlated with the deviation 
from the LV relation at the central region. In this sim¬ 
ulation, the kinetic energy is not sufficiently gained to 
attain the LV relation, which reflects both the tempera¬ 
ture falling off and the steepness of the phase space den¬ 
sity which is derived from dividing the density with the 
falling temperature. 


V. CONCLUSIONS AND DISCUSSIONS 

We focused on the local virial (LV) relation and investi¬ 
gated the role of the LV relation in a spherically symmet¬ 
ric steady state system. Assuming the constant LV ratio 
and solving Jeans equation for spherically steady state 
system, we studied the full solution space of the prob¬ 


lem under constant anisotropy parameter and obtained 
some meaningful solutions. Especially, if the LV relation 
is satisfied, the solution of the power law behavior of the 
density profile in both the asymptotic regions r —> 0 and 
oo exists and is commonly observed in many numerical 
simulations. In this sense, the LV relation plays a critical 
role in solutions of Jeans equation for spherically steady 
state system. 

In the previous section, we compared our cold col¬ 
lapse simulations with the critical solution with con¬ 
stant anisotropy parameter (3. We got the results that 
the simulations agree with the critical solutions quite 
well except for the initially high density simulation with 
(b t , a) = (0,2). In that simulation, the LV ratio b be¬ 
comes smaller than one around the center, where the tem¬ 
perature falls off and the phase-space density becomes 
steeper than the mass density. 

These results remind us of the cosmological simulations 
for which both mass density and the phase space density 
have the universal characters. In the cosmological simu¬ 
lation, the phase space density behaves as p/a 3 ~ r _ “, 
where a ~ 1.87 0,113' which is larger than the exponent 
of the mass density in the central part and the temper¬ 
ature falls off in the central region. Hence we speculate 
that both of the initial high density simulations and cos¬ 
mological simulations have the common characters for 
the phase space distribution. This may be because in 
the initially high density simulation, the matters in the 
central part collapse to make a core in the early stage 
and the matters in the outer part gradually falls into the 
core in the later stage, which is similar to the case with 
the secondary in-fall in the cosmological simulation. 

However, we have to be careful to conclude that these 
results mean the inapplicability of the LV relation for 
these cases, because we compared these cases only un¬ 
der the condition (3 is constant. It is not obvious that 
the same kind of asymptotic behavior of the critical so¬ 
lutions obtained under the LV relation are sustained for 
any function form of /3(r). We will discuss these points 
in our coming paper. 


APPENDIX A: A SOLUTION FOR a = 1 CASE 

1. (3=1 case 

In the case for (/3, b) = (1,1), the equation © becomes 

i = ( 7 - 7 +X 7 - 7-) = «> (Al) 


where 


7 + = 


5 + VI — 4k 


7 _ = 5-vT^ m 


Integrating the Eq. ED- we obtain 


7 + + 7_a; 7 + 7 ~ 

l + z 7 +- 7 - 


(A3) 











where we choose the integral constant as satisfying the 
condition 7(1) = (7+ +7_)/2. At large x, 7 approaches 
7_ and the total mass is infinite because 7_ < 5/2. From 
Eqs.0, 0, and I A: 11 ) . we have 

p oc x _7+ (l + x 7+-7 -) , (A4) 

( -1 <y_)_—_ \ 

-h —-) , (A5) 

3-7+ 3 - 7_ J 

</> oc x 2_7+ (l + a; 7+_7_ ) , (A6) 

er 1 2 oc -x 2_7+ (l + ir 7+-7 ^) , (A7) 

p/a 3 oc 2 3/2 x 7+/2 “ 3 (1 + x 7 +" 7 -)" 1/2 , (A8) 

where M := Jg 47 xu 2 pdu. 

2. /3 ^ 1 case 

In this case, the equation becomes 

7'- (7 - 7a)(7 - 7b) = kx 2(1_/3) , (A9) 

where 

7a = 2/3, 7 b = 2/3 + 1. (A10) 

The behavior of the solutions in Eq. shows in Figca 
In the limit r —> 0, 7 approaches "f a or 7f, or —00 (inner 
hole). On the other hand, in the limit r —> 00, all so¬ 
lutions have an outer truncation at r = 00. Physically 
meaningful solution is represented by a solid line in FigH5l 
where 7 approaches 7 a at r —> 0. 



FIG. 15: The same as Fig0J but for case a = 1 (/? ^ 1). The 
parameters are set to (/?, 6) = (3/4, 3/2). 


APPENDIX B: THE CONNECTION OF 
RELATIVE POTENTIALS 

Here we connect two critical solutions with different 
values of constant /3, which leads to the critical solution 
with a step function form of /3 (r). Here we use the unit 
G = r c = (f> c = 1, where r c is the radius of the connected 
position and </> c is the value of the relative potential at 
r = r c . In this unit the relative potentials <f> 1 in the inner 
region and </>2 in the outer region are described as 


II 

( l + ci 1 N 

l/si 

(Bl) 

\r Sl + c® 1 , 

) 5 

/>2 = 

( 1 + 4 2 ' 

1/S2 

(B2) 

U S2 +C 2 2 , 

j ’ 


respectively. Each potential includes a free parameter ci 
or C2 other than the parameter Si and S2 related to the 
anisotropy of each region. From the continuous condition 
of the first derivative of </>+ and </>2 against r at r = 1, C2 
is described as 


c 2 =c® 1/S2 . (B3) 

Hence connected solution depends on the one-parameter 
c = C\ for given s 1 and s 2 . 

This parameter c determines the mass fraction of the 
inner solution </> 1 against total mass as 

TF- = (l + c si r 1-1/s2 . (B4) 

Mtot 

Hence we can adjust the mass fraction of the inner region 
by changing the parameter c. 


ACKNOWLEDGMENTS 

The authors would like to thank Professor Kei-ichi 
Maeda for the extensive discussions. All numerical sim¬ 
ulations were carried out on GRAPE system at ADAC 
(the Astronomical Data Analysis Center) of the National 
Astronomical Observatory, Japan. 


[1] M. Henon, Ann. Astrophys. 27, 83 (1964). 

[2] T. S. van Albada, Mon. Not. R. Astron. Soc. 201, 939 
(1982). 

[3] J. Binney and S. Tremaine, Galactic Dynamics (Prince¬ 

ton University Press, Princeton, 1987). 


[4] O. Iguchi, Y. Sota, T. Tatekawa, A. Nakamichi, and 
M. Morikawa, Phys. Rev. E 71, 016102 (2005). 

[5] Y. Sota, O. Iguchi, M. Morikawa, and A. Nakamichi, 
astro-ph/0403411 submitted to Phys.Rev.E, (2005). 

[6] N. W. Evans and J. An, Mon. Not. R. Astron. Soc. 360, 









9 


492 (2005). 

[7] A. Kawai, T. Fukushige, J. Makino, and M. Taiji, 
Publ. Astron. Soc. Japan 52, 659 (2000). 

[8] W. Dehnen and D. E. McLaughlin, Mon. Not. R. As¬ 
tron. Soc. 363, 1057 (2005). 

[9] Y. Sota, O. Iguchi, M. Morikawa, and A. Nakamichi, Pro¬ 


ceedings of CN-Kyoto International Workshop on Com¬ 
plexity and Nonextensivity, (2005). 

[10] U. -I. K. Veltmann, Astron. Zh. 56, 976 (1979). 

[11] J. E. Taylor and J. F. Navarro, ApJ. 563, 483 (2001). 



